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By considering adiabatic contraction of the dark matter (DM) during the star formation, we 
estimate the amount of DM trapped in stars at their birth. If the DM consists partly of primordial 
black holes (PBHs), they will be trapped together with the rest of the DM and will be finally 
inherited by a star compact remnant — a white dwarf (WD) or a neutron star (NS), which they will 
destroy in a short time. Observations of WDs and NSs thus impose constraints on the abundance of 
PBH. We show that the best constraints come from WDs and NSs in globular clusters which exclude 
the DM consisting entirely of PBH in the mass range 10 16 g — 3 x 10 22 g, the strongest constraint on 
the fraction ^pbh /^dm J$ 1CP 2 being in the range of PBH masses 10 17 g- 10 18 g. 



PACS numbers: 

I. INTRODUCTION 

Various observational evidences point at the existence 
of a new matter component in the Universe, the dark 
matter (DM) (for a recent review see, e.g., [HH]). Obser- 
vations of the cosmic microwave background imply that 
DM comprises ~ 23% of the total energy budget of the 
Universe, thus dominating in the matter sector, where 
baryonic component sums up to only 4% [3]. However, 
the nature of DM remains unknown and masses of pos- 
sible candidates range over many orders of magnitude 
from a fraction of eV to many solar masses. Although 
most popular candidates are new stable particles, other 
possibilities are not excluded. 

In the Early Universe density perturbations with high 
initial amplitude could collapse forming black holes [3]. 
If some of these black holes survive until now they could 
constitute (at least) a fraction of DM. Properties of these 
primordial black holes (PBHs) make them a suitable DM 
candidate: they are nonrelativistic and have subatomic 
size r ~ 10~ 8 cm (tobh/10 2O <7) which makes them effec- 
tively collisionlcss. Unlike most of other DM candidates, 
PBHs do not require the existence of new particle species. 

The initial mass function of PBHs is flat in the case of 
a flat power spectrum of primordial density fluctuations. 
However, models with strongly non-flat mass function of 
PBHs can be constructed, see, e.g., Refs. [B]. The 
constraints at different masses, therefore, should be con- 
sidered independently. 

The PBHs with masses tobh < 5 x 10 14 g have lifetimes 
due to Hawking evaporation [7 shorter than the present 
age of the Universe. Such PBHs thus cannot contribute 
to the DM. 

PBHs with slightly larger masses emit 7— rays with 
energies around ~ 100 MeV [8]. Observations of the 
extragalactic gamma-ray background with the Energetic 



Gamma Ray Experiment Telescope (EGRET) [§] set an 
upper limit on the cosmological density Opbh of such 
PBHs as a function of their mass, e.g. ^pbh < 10~ 9 for 
tobh = 10 15 g [TO]. These observations show that PBHs 
with masses Tobh < 10 g can not constitute more than 
1% of DM. However, the constraints coming from the pro- 
cess of Hawking evaporation disappear for PBH masses 
larger than tobh > 7 x 10 16 g. 

The PBHs in the mass range tobh ;$ 10 19 — 10 20 g 
can be constrained with the so-called "femto-lensing" of 
the gamma-ray bursts |llj . Present day observations of 
gamma-ray bursts constrain the mass fraction of PBHs 
in the narrow mass range around tobh ~ 10 18 g at sev- 
eral percent level [12] . The abundance of more massive 
PBHs can be constrained from microlensing surveys. The 
EROS microlensing survey sets an upper limit of 8% on 
the fraction of the Galactic halo mass in the form of PBHs 
with masses in the range of 10 26 g < tobh < 3 x 10 34 g 
[13]. At even higher mass scales, 10 33 g < tobh < 10 40 g, 
the analysis of the Cosmic Microwave Background can 
be used to constrain PBHs at the level of 10 -7 [T4"] . 

The range of PBH masses from roughly 10 17 g to 10 26 g 
remains essentially unconstrained, apart from the above- 
mentioned narrow region around 10 18 g. The aim of this 
paper is to constrain PBHs as the DM candidates in this 
still allowed mass range. To this end, we investigate the 
effect of PBHs on the evolution of compact stars - neu- 
tron stars (NSs) and white dwarfs (WDs). The main 
idea is as follows. PBHs may be captured by a star in 
the process of its formation. This leads to no observa- 
tional consequences until the evolution of the star reaches 
the stage of a neutron star or a white dwarf. Then the 
accretion onto a PBH becomes sufficiently fast to destroy 
the compact star in a short time [T51 - H7] . The region of 
PBH masses and abundances where this happens with 
large probability is thus excluded by observations of the 



2 



existing neutron stars and white dwarfs. 

The paper is organized as follows. In Sect. II we dis- 
cuss the gravitational capture of DM during the process 
of star formation. In Section III we derive the constraints 
on the fraction of PBHs in the DM from the existence of 
WDs and NSs. In Section IV we summarize the results 
and present our conclusions. Throughout the paper, we 
use the units h = c = 1 . 



II. CAPTURE OF DARK MATTER DURING 
STAR FORMATION 



In this section we study the capture of dark matter 
during the star formation process, neglecting all the DM 
interactions except the gravitational one. The purpose 
is to estimate the total amount of DM captured inside 
newly formed stars. 



A. Star formation stages 

Star formation occurs mainly in giant molecular clouds 
(GMCs). GMCs are dense regions of the interstellar 
medium composed primarily of molecular hydrogen (H 2 ) 
with typical mass M ~ 3 x 10 5 Mq and average density 
p ~ 550 GeV cm -3 , which would imply a radius of 17pc 
in the case of a spherical shape. A GMC is usually com- 
posed of smaller overdense sub-clouds, i.e., clumps. In 
gravitationally bound cores inside the clumps, individual 
stars are formed. 

The formation of stars involves different stages. The 
first one corresponds to the fragmentation of a GMC into 
gravitationally bound regions that are initially supported 
against gravity by a combination of rotation, magnetic 
and turbulent pressures [T51 [TH] - At some point, as the 
cloud core loses its magnetic and turbulent support by 
still poorly understood mechanisms like ambipolar diffu- 
sion [IS], the growing central concentration becomes un- 
stable to the gravitational collapse. At the initial stage 
of the collapse, the cloud has a uniform temperature, is 
rotating slowly and has an almost flat density profile in 
the central part. At the end of this phase an opaque pro- 
tostcllar object in hydrostatic equilibrium is formed at 
the center [20H22]. 

When the protostar is formed, it accretes from the sur- 
rounding disk increasing its temperature. When the cen- 
tral object has accumulated most of its main-sequence 
mass and the surrounding disk disperses, it is considered 
a pre-main-sequence star. The main energy source for 
such an object is the gravitational contraction, contrary 
to nuclear fusion for main sequence stars. Therefore it is 
evolving on a Kelvin-Helmholtz timescale GM% /{R+L*), 
where M^^R* and are the mass, the radius and the 
luminosity of the pre-main-sequence star, respectively. 
This time is longer than the free fall time (R^/GM*) 1 ' 2 . 



B. Adiabatic contraction 

The main mechanism of the capture of DM by stars at 
the time of formation is the adiabatic contraction. Con- 
sider first this mechanism in general terms. 

In this paper we will be interested in systems that are 
dominated by baryons. In this case the adiabatic contrac- 
tion is easy to understand. When baryons contract los- 
ing energy by non-gravitational mechanisms, their time- 
dependent gravitational potential pulls the DM particles 
along. The DM distribution thus develops a peak cen- 
tered at the baryon distribution. 

If the change of the baryonic gravitational potential is 
slow, that is, if the characteristic time of the baryonic 
contraction is much larger than the free fall time iff — 
(Gpo)~ 1/2 , where Po is the baryonic density of a cloud, 
the DM distribution is determined by the (approximate) 
conservation of the adiabatic invariant 



pdq = ET 



(1) 



where p and q are the phase space coordinates of a DM 
particle of energy E and orbital period T. Moreover, the 
angular momentum is conserved for each DM particle as 
long as the potential is central. From these conserved 
quantities, a relation between the initial orbital radius 
and the final one can be derived [2"3Tf2T)] . 

Regardless of whether the contraction of DM is adia- 
batic or not, the phase-space density of DM has to be 
preserved all along the contraction process, as dictated 
by Liouville's theorem. For an initial Maxwellian velocity 
distribution of DM with the dispersion v, the maximum 
phase space density is at zero velocity and equals [27] 
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where /?dm is the space density of DM and todm is the 
DM mass. The effect of the adiabatic contraction is to 
fill all the allowed phase-space with the density close to 
the maximum value. 

In the case of circular orbits the conservation of the an- 
gular momentum and the adiabatic invariant (JlJ implies 
the conservation of the quantity rM(r), where M(r) is 
the mass within the radius r. Suppose a baryonic cloud, 
which was initially a uniform sphere of radius R, con- 
tracts to a compact object of a negligible size. Assuming 
that the DM particles move on circular orbits, the initial 
uniform DM density pdm is modified as follows [28] : 



Pdm(V) = -Pdm 



1_ fR\ 
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J 



(3) 



provided the adiabatic approximation holds. 

In realistic cases the orbits of DM particles are not 
circular. The question thus arises whether eq. (|3| is a 
good approximation in realistic situations. An exact cal- 
culation has been performed in another limiting case of 
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FIG. 1: DM density profile obtained from the simulation 
after the adiabatic formation of a star. The inner part of the 
profile has the slope close to —3/2 as expected from the Li- 
ouville's theorem for the uniform initial velocity distribution. 



density profile ([3]), even for a relatively rapid change of 
the external potential, t c — 3 . In the case of random 
initial velocities, however, the inner profile was found to 
have a slope close to —3/2, 



M 1 ~ (R 



3/2 



(4) 



as represented in Fig. [T] These results are in agreement 
with Liouville's theorem, since the final DM velocity goes 
as v(r) oc r -1 / 2 . Since random initial velocities appear to 
be a better approximation to realistic initial conditions 
than the circular ones, and because the profile Q gives 
more conservative estimates, in the rest of this paper we 
use the profile Q. 



C. DM bound to a baryonic cloud 



purely radial orbits [551 EH]- The result were found to 
be roughly compatible with the case of circular orbits. 
However, this is not a realistic case either. Another ques- 
tion of practical importance is the domain of validity of 
the adiabatic approximation. Formally, it requires the 
time of collapse t c to be much longer that the free-fall 
time tff. In Ref. [3T] high-resolution numerical simula- 
tions have been performed and it has been shown that 
the adiabatic contraction may remain a good approxima- 
tion, depending on the potential, even when t c ~ iff . 

To clarify these issues we have performed the following 
numerical simulation. For the baryonic distribution that 
is responsible for the time-varying external gravitational 
potential we took the sum of a uniform spherical cloud 
and a point mass in the center. The point mass was zero 
at the initial moment of time and then increased linearly 
with time, while the mass of the spherical part, always 
uniform in density, decreased in such a way that its sum 
with the point mass remained constant. The time t c over 
which all the mass was transferred from the cloud to the 
central object was treated as a free parameter. 

The DM particles were injected at t = with an initial 
uniform distribution in position and velocity. The initial 
density profile was taken to be constant over the volume 
of the cloud to mimic the physical properties of prestellar 
cores. The particles were injected one by one, which cor- 
responds to neglecting the DM contribution to the grav- 
itational potential. Those particles which have positive 
total energy at t = (and thus are not gravitationally 
bound to the system) were discarded; the remaining ones 
were evolved numerically in the time-dependent gravita- 
tional potential. At a random time t > t c the positions of 
these particles were sampled in order to reconstruct, after 
many simulations, the final density profile. As a consis- 
tency test, we have also performed an identical simulation 
with initial velocities of DM particles generated in such 
a way that the particles move on circular orbits. 

In the case of circular orbits we have reproduced the 



As is clear from the above discussion, only DM grav- 
itationally bound to a baryonic cloud is subject to the 
adiabatic contraction when the cloud collapses. There- 
fore, to set the initial conditions for the adiabatic con- 
traction we need to estimate the amount of DM that is 
gravitationally bound to a cloud. 

We will assume that the matter density is dominated 
by baryons, as is the case in the star forming regions. 
When the overdensity of baryons is formed, some amount 
of DM ends up gravitationally bound to the baryonic 
cloud. Consider the formation of a spherical cloud of ra- 
dius Rq and baryonic density p . Our goal is to estimate 
the density of DM bound to the cloud, pdm. bound, given 
the mean density of DM, /?dm- We will assume that orig- 
inally the DM particles have the Maxwellian distribution 
in velocities with the dispersion v, 
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where Hum = Pdm/w, m being the mass of the DM 
particle. We will see that in the cases of interest the 
velocities of bound particles are much smaller than v, and 
thus the precise shape of the distribution is not essential. 

After the formation of a baryonic cloud, the gravita- 
tional potential felt by DM particles becomes of order 

4, ~ O = 2nGp R 2 . 

Those particles with kinetic energies smaller than cf>o 
(equivalently, velocities v < vq = y/2(po) become grav- 
itationally bound. Their number density is obtained by 
integrating eq. (|5| up to v = vo. Multiplying by the DM 
mass, one has 
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M*/M 


po, GeVcm~ 3 


Ro, AU 


1 


10 6 


4300 


2 


1.8 x 10 6 


4450 


3 


2.4 x 10 6 


4620 


4 


3.1 x 10 6 


4710 


5 


3.6 x 10 6 


4780 


6 


4.2 x 10 6 


4840 


7 


4.8 X ID 


4880 


8 


5.3 x 10 6 


4930 


10 


6.4 x 10 6 


5000 


12 


7.4 x 10 6 


5060 


15 


8.8 x 10 6 


5130 



TABLE I: The parameters of prestellar cores used in the 
estimates. 



where we have assumed vq <C v and thus set the expo- 
nential to 1. 



D. Globular clusters 

Globular clusters (GCs) are gravitationally bound sys- 
tems consisting of 10 4 to 10 7 stars with average diam- 
eters ranging from 20 to 100 pc. There are about 100 
GCs known in our Galaxy. A typical GC has a baryonic 
mass of (a few) x 10 5 A/ Q and a core radius of 1 — 2pc. 
The age of GCs is about 8 to 13.5 Gyr [32], and as such 
they are the oldest surviving stellar subsystems in the 
galaxy, made up of the population II stars, white dwarfs, 
neutron stars and black holes. 

There are two classes of scenarios for GC formation. 
According to the primordial, or "DM-dominated" one, 
GCs were formed by the infall of baryonic matter into the 
gravitational wells of the DM density peaks at redshifts 
z > 10 [33 38J . The second is the "baryon-dominated" 
scenario, according to which GCs were formed in predom- 
inantly baryonic processes like major mergers, hydrody- 
namical shocks and so on mostly in proto-galaxies that 
later assembled into the Milky Way [3^142) . Moreover, 
both mechanisms could be at work because the observed 
distribution of metallicity in GCs is clearly bimodal, so 
that metal-poor GC could be of cosmological origin, while 
metal-rich GCs could be formed in the course of mergers 
[33] . Although there is no evidence of DM presence in the 
GCs now [33], it was shown that it could be present at 
the formation time and subsequently tidally stripped due 
to interactions with the host galaxy [35] . We will assume 
in what follows that at least some of the GCs resided in 
DM minihaloes in the past, and concentrate on those. 

The DM density in the central regions of GCs has been 
estimated in Ref. |46j by making use of the formalism 
developed in Refs. [33 133- The conclusion was that the 
present-day DM density close to the core of a GC is of or- 
der p DM <~ 2 x 10 3 GeV cm~ 3 , the estimate being rather 
insensitive to the original halo mass. This result is in con- 



cordance with the N-body simulations [55J [3SJ suggesting 
that the inner part of DM halos survive successive tidal 
interactions with the host galaxy. As has been stressed in 
Ref. [3S], the number cited includes the effect of dynam- 
ical heating of DM by the stars comprising the cluster, 
which reduces the DM density. In our estimates this ef- 
fect is irrelevant since we are interested in the evolution 
stage prior to the star formation. With no heating, the 
DM density in the core would be pdm ~ 10 4 GeV cm~ 3 , 
which we adopt in what follows. 

Another important parameter is the value of the DM 
velocity dispersion in GCs. As stars in the GC are col- 
lisionless and behave essentially as DM particles, this 
parameter can be extracted directly from observations. 
Although there is quite a bit of scatter, typical observed 
GCs have the velocity dispersion around v — 7 kms -1 

m 

Let us now turn to the prestellar cores. Their typical 
parameters are known from observations carried out with 
the SCUBA instrument [48]. The dataset is well fitted by 
the Bonnor-Ebert profile [3S] which, for our purposes, can 
be approximated by the flat core of radius Rq containing 
the baryonic mass M* of the future star. 

Two cases will be of interest in what follows: stars with 
masses 1M© < M* < 7M Q which are typical progenitors 
of a WD, and supermassive stars of masses, M* > 8M 
progenitors of NSs. In all cases the gravitational poten- 
tial of the prestellar core is (much) smaller than that of 
the GMC, so one may use again eq. ^ to estimate the 
density of DM gravitationally bound to the core. We use 
the parameters of prestellar cores that are listed in Ta- 
ble [I] As has been already mentioned, the formation of 
the prestellar cores relies on the non-gravitational energy 
loss mechanisms and thus is expected to be slower than 
the free fall. 

Making use of eq. Q , one obtains the total DM mass 
contained in a star formed within the GCs as listed in Ta- 
bic [Tl] These values were calculated with the DM density 
and velocity dispersion given above; for different values 
of these parameters the mass of the bound DM should 
be rescaled by the factor 

( - \ 3 ( ^ DM ^ ( 8 ) 
V7km/ S y \ v 10 4 GeV/cm 3 y 

which may be different for different GCs. 

III. CONSTRAINTS ON PRIMORDIAL BLACK 
HOLES 

So far the discussion has been general and did not de- 
pend on the DM nature. Consider now specifically the 
case of primordial black holes. The PBHs that end up in- 
side a star when the latter is formed start accreting and 
gravitationally pulling on the surrounding matter, lose 
their momentum and gradually sink to the center. The 
sinking process is slow, so that the characteristic time 
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M*/M Q 


Ppsc, GeVcm 3 


Afbound, g 


1 


2 x 10 1 


4.4 


x 10 19 


2 


5.2 x 10 1 


2.5 


x 10 20 


3 


9.2 x 10 1 


7.2 


x 10 20 


4 


1.4 x 10 2 


1.5 


x 10 21 


5 


1.9 x 10 2 


2.6 


x 10 21 


6 


2.4 x 10 2 


4.2 


x 10 21 


7 


3 x 10 2 


6.2 


x 10 21 


8 


3.6 x 10 2 


8.7 


x 10 21 


10 


5 x 10 2 


1.6 


x 10 22 


12 


6.4 x 10 2 


2.4 


x 10 22 


15 


8.7 x 10 2 


4.3 


x 10 22 



TABLE II: Density of DM bound to the prestellar core, 
Ppsc, and the total mass Mb OU nd of DM contained in a star 
right after its formation in a GC with the central DM density 
Pdm ~ 10 4 GeV cm and velocity dispersion v = 7 kms 
for different star masses. 




' le+15 le+16 le+17 le+18 le+19 

BH mass, g 

FIG. 2: The dependence of the size r c of the collection region 
(the region from which the PBHs captured by the star at its 
formation have enough time to sink to within the radius of 
the future compact remnant, WD or NS) on the PBH mass, 
corresponding to the case of WD for M„ = Mq . 



may exceed the age of the star. Because of their slow 
accretion and small number, the presence of BH has no 
observable effects on the star evolution at this stage. 

When a star polluted by PBHs evolves into a compact 
object (WD or NS), some of the BHs get inside the com- 
pact remnant. Because of a much higher density, the 
accretion is now more efficient and the PBHs, if present 
inside the remnant, rapidly consume the latter. The ob- 
servation of WDs and NSs thus implies constraints on 
the abundance of PBHs, which has to be such that the 
probability to get a PBH inside NS or WD is much less 
than 1. 

To quantify this statement we calculate the number 
-^bh of BHs that would sink down to the future radius 
Tf of the compact remnant by the end of the star evolu- 
tion and, thus, would end up inside the star remnant. If 
-/Vbh < 1 no constraints arise. If TVbh > 1, the maximum 
allowed fraction of BHs in the total amount of DM is 



^pbh < 1 



Q 



DM 



BH 



(9) 



Thus, in the range of PBH masses where -/Vbh > 1, PBHs 
cannot constitute all of the DM. 

The PBHs that are eventually trapped by the compact 
remnant initially occupy some spherical volume of radius 
r c which we call the "collection region" . Knowing r c as a 
function of the PBH mass 771bh and the DM distribution 
inside the star at the time of formation allows one to 
calculate ./Vbh as follows, 



TVbh = AfcitW/mBH, 



(10) 



where Mtjm (r) is the DM mass contained in the radius r 
at the time of the star formation. 

The sinking of the PBH inside the star has been con- 
sidered in Ref. [30] ■ The dynamical friction force per unit 
PBH mass is given by the Eq. (16) of Ref. [3D]. Multi- 
plied by the PBH velocity, this gives the PBH energy 



loss rate dE/dt. On the other hand, assuming circular 
orbits, dE/dt can be expressed in terms of the change of 
the orbit radius. Equating the two gives a closed first 
order differential equation for the orbit radius as a func- 
tion of time, r(t). We derive and solve the corresponding 
equation numerically in Appendix [XJ assuming the star 
model with the polytrope index n — 3. 

Having found the dependence r(t) for a given BH mass, 
we fix the final radius rt = r(i*) to be the size of the 
compact object (NS or WD). Here t* is the lifetime of 
the star. We then determine the collection radius as 
r c = r(0). (In practice, it is more convenient to run 
the evolution equations backwards in time starting from 
r = r / .) The dependence of r c on the BH mass is shown 
in Fig. [2] for A/* = M Q and r; = r WD = 10 4 km. 

At small BH masses the dynamical friction is inefficient 
and the collection radius r c is not very different from 
rt. As the BH mass gets larger the friction becomes 
more efficient, so that r c grows and eventually becomes 
close to the star radius. The transition is quite rapid; 
the value of tobh = 77Hrans at which it occurs can be 
understood analytically from the behavior of eq. ( A8 1 , 



see Appendix[A]for details. It corresponds to the smallest 
ttibh for which the collection radius is still close to the 
star size (i.e., the lifetime of the star is still sufficient for 
a BH to sink to r = rt starting near the surface). By an 
order of magnitude, it is given by 



"Itrans ~ 4 X 10 



10 Gyr 



1/2 



R* 

R Q 



3/2 



(11) 



Note that there is no dependence on 77 because the final 
stages of the BH sinking are exponential, and these are 
the (longer) initial stages that set the overall time scale. 

The number of BHs inside the collection region -/Vbh 
can be found from the total DM mass trapped by the 
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FIG. 3: Constraints on the fraction S^pbh/^dm • Purple 
shaded region is excluded by observations of WDs and NSs in 
the centers of globular clusters. Thin curves show the exclu- 
sions from different star masses. 



star (see Table [nj) and the DM density profile inside the 
star. The DM profile inside the star after the adiabatic 
contraction is determined by eqs. ([3| and Q and the 
baryonic density, for which we assume the density profile 
of the polytrope n = 3 model. From cq. Q , the DM and 
baryonic masses are related as follows, 



MdmO") = A'f bound 



( M(r)r c 



1/2 



(12) 



where M(r) and A/dmM are the baryonic and DM mass 
within the radius r, respectively. The number of BH 
within r c is then given by eq. (10 1. 



The resulting constraints on the fraction of PBHs in 
the total amount of DM are shown in Fig. [3] Purple 
shading shows the region excluded by the observations 
of WDs and NSs in the globular clusters. Thin curves 
show the exclusion regions resulting from different star 
masses. One can see that the constraints from WDs and 
NSs complement each other and together cover the range 
of masses from 10 16 g to 3 x 10 22 g. 

The shape of the excluded regions is similar in all cases 
shown in Fig. [3] It can be understood from the mass 
dependence of the collection radius r c , Fig. [2] as follows. 
At large PBH masses the size of the collection region is 
close to the star size, so that Mdm(^c) — Abound and 
the maximum PBH fraction Opbh/^dm scales like ttibh, 
i.e., the constraints improve at smaller masses. However, 
at some point around tobh ~ Wtrans the collection radius 
r c decreases (cf. Fig. [2]) and the constraints get relaxed 
very rapidly. 



IV. CONCLUSIONS 



would pollute the newly-formed stars and, after sinking 
to the center, would end up in the compact remnant re- 
sulting from the star evolution (WD or NS). Once inside 
the remnant, PBHs would rapidly destroy it by accretion. 
Mere observations of WDs and NSs, therefore, impose 
constraints on the abundance of PBHs. 

We have found that the most stringent constraints 
come from observations of WDs and NSs in globular 
clusters. WDs and NSs are sensitive to the mass ranges 
10 16 g < m BH < 10 21 g and 10 21 g < to B h < 3 x 10 22 g, 
respectively, thus complementing each other. Every- 
where in this mass range the PBHs are excluded as com- 
prising all of the DM. The best constraint on the PBH 
fraction Hpbh/^dm J$ 10 -2 was found for ttibh in the 
range 10 17 g — 10 18 g. 

The constraints derived from the globular clusters are 
based on the assumption that at least some of those were 
formed in a primordial DM-dominated environment. As 
a word of warning, it should be noted that this issue is 
still debated in the literature. For instance, observations 
of a low-metallicity cluster NGC 2419 [49] [50] seem to 
indicate that its mass-to-light ratio is in a good agree- 
ment with what is expected for a pure baryonic system. 
However, NGC 2419 has a number of extreme properties 
[5TH53] that make the globular cluster nature of this ob- 
ject questionable. In addition, high-resolutions N-body 
simulations [35] SH] indicate that the mass-to-light ratio 
may not be sensitive to the presence of the DM compo- 
nent in GCs. 

In order to derive the constraints on the PBH abun- 
dance we have investigated the baryonic contraction of 
the DM during the star formation process. In particular, 
we have calculated numerically the resulting DM profile 
and found the slope close to —3/2. We also estimated 
the total amount of DM that is trapped inside the star 
at the time of its formation. This part of our results is 
not specific to any particular form of the DM. 
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We have derived the constraints on the abundance of 
PBH from observations of the existing WDs and NSs. 
The origin of these constraints is as follows. If PBHs were 
present at the time of star formation, i.e., at z < 10, they 



Appendix A: Sinking of BH in the star 

Here we consider the energy loss by a BH that is or- 
biting inside a star gradually sinking to the center. The 
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star is assumed to have an r-dependent density p(r) and 
temperature T(r), and the mass M(r) enclosed inside the 
radius r. Our purpose is to derive the equation that gov- 
erns the evolution of the BH orbit due to the dynamical 
friction, assuming the orbit is circular and changes slowly. 
The question of dynamical friction was considered in a 
general context , e.g., in Refs.[54"l 155]. 

The BH moving through a star experiences a dynami- 
cal friction force [5S] that can be written as follows, 



= -7(w)v, 



(Al) 



where 



-y(u) = 47rG 2 p(r)TO B Hln(A) 



F(X) 



F(X) = erf (A) - 2X exp(-A 2 )/ % A : , (A2) 



X = v/(V2a), 
p{r) is the density of the baryonic gas comprising the star, 



a is the velocity dispersion of the particles a — yjT/m 
with temperature T and mean molecular weight m 
(m ~ 1.6 GeV for a main sequence star), and ln(A) ~ 
ln(M*/mBH) — 30 is the Coulomb logarithm [SSJ. Multi- 
plying Eq. by v gives the total BH energy loss rate 
per unit mass, dE/dt = —j(v)v 2 . 
Making use of the relation 



= GM(r)/r, 



(A3) 



the same energy loss rate can be expressed through the 
change of the radius r of the BH orbit, 



dE 
~dt 



d 
dt 



1 



U(r) 



dr v 2 ( 4irr 



dt 2r { M(r) 



p(r) + 1 



where U(r) is the gravitational potential. Equating the 
two quantities and simplifying by v 2 one gets 



dr 1 ( Anr 3 
dt2r \M(r) 



p{r) + 1 > = -j(v). 



(A5) 



By virtue of eq. ( A3 ), this is a closed differential equation 
for the BH orbit radius r(t) as a function of time. Note 
that this equation is not equivalent to eq. (24) of Ref. [30] 
because the contribution of the gravitational potential 
(the term U{r) in eq. (A4)) has been missed there. 

Let us rewrite this equation in the form convenient for 
the numerical solution. Define the dimensionless quanti- 
ties 



x 

T 

p(r) 
M(r) 
T(r) 



r/R*, 
t/t , 

p(r)/p(0), 
M(r)/M„ 
r(r)/T(0), 



where 



M, 



3/2 



27iVG / 9(0)m B H#* /2 hi A 



4.2 x 10 3 yr 



mBH_ 

lQ22g 



1/2 



Rq 



.3/2 



(A6) 



(A7) 



The profiles p(r) and T{r) are determined by the star 
model. The normalization parameters /o(0), T(0), M* 
and are not independent. They obey the following 
two relations, 



GM»m 
R,T(0) 

p(Q)R 3 * 



= 1.17 



12.9, 



where, as before, m is the mean molecular weight. The 
scaling in eq. (A7| takes into account these relations. 



In terms of the dimensionless quantities, Eq. (A5| be- 
comes 



dx 



p(x) 



f(x)M 3 / 2 (x 
Here we have introduced the function 



F(X). 



(A8) 



m = \ 



4tt 



r 3 p(r) 
M(r) 



1 



4 1 M(x) 



which varies between 1 in the star center x = and 1/4 
at the star surface x = 1. The variable X is in turn a 
function of x which can be expressed as follows, 



X- Y /J = 0.765 



2rT(r) 




while the function F(X) is defined in Eq. ( A2 1. At small 
values of X this function behaves as 4X 3 / (3i/7r); at large 
X it asymptotes to 1. 



A useful analytical insight into behavior of the eq. ( A8 ) 
can be obtained by considering two limiting cases. "At 
small values of x such that the parameters of the star 
can still be approximated by their core values one has 

X ~ 5.62x. 



At small x such that X 1 eq. ( A8 ) becomes 
dx 



dr 



= -0.337a;, 



whose solution is x(t) = cxp(— 0.337t/t ) with i given 
by eq. ( |A6[ ). In this regime, valid for the final approach 
by a sinking BH of the radius rf (the radius of a future 
compact object), the characteristic time scale is 

(At) 2 ~3\n(r /r f )t Q , 



8 



where ro is some initial radius. 

At moderately small x such that X > 1 eq. ( A6 ) takes 
the form 



dx 



1 



397 a; 2 ' 

which gives the evolution time (At)! from x\ to x^ 



where we have set x\ ~ 1 and neglected x\. We see 
that this first stage is typically longer than the second, 
(Ai)i > (At)2- Equating (At)i to the lifetime of the star 
and making use of eq. ( A6 ) leads to the estimate (111. 



(At)i = 132 f (a;? - x\) ~ 10 2 x i 
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